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Abstract 


This report describes a new method for determination of the geopotential, or the 
equivalent geoid. It is based on Satellite-to-Satellite Tracking (SST) of two co-orbiting low 
earth satellites separated by a few hundred kilometers. The analysis is aimed at the GRACE 
Mission, though it is generally applicable to any SST data. It is proposed that the SST be viewed 
as a mapping mission. That is, the result will be maps of the geoid or gravity, as contrasted with 
determination of spherical harmonics or Fourier coefficients. A method has been developed, 
based on Geophysical Inverse Theory (GIT), that can provide maps at a prescribed (desired) 
resolution and the corresponding error map from the SST data. This computation can be done 
area by area avoiding simultaneous recovery of all the geopotential information. The necessary 
elements of potential theory, celestial mechanics, and Geophysical Inverse Theory are described, 
a computation architecture is described, and the results of several simulations presented. 
Centimeter accuracy geoids with 50 to 1 00 km resolution can be recovered with a 30 to 60 day 
mission. 
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I 


Introduction 


This analysis is focused on the Gravity Recovery and Climate Experiment (GRACE) 
mission, to map the earth’s gravity field, or geopotential. The GRACE mission concept was 
based on extensive experience with dynamic satellite geodesy, using well established methods 
and with the confidence based on very good results. Though there is no doubt that the existing 
methods will satisfy the mission requirements, we propose here an alternative methodology for 
analyzing the mission data that can provide a different type of result, some analysis advantages, 
and possibly lower analysis cost. The basic objective is to determine regional geoid maps, 
directly from the satellite-to-satellite (SST) tracking data, without recourse to global spherical 
harmonic solutions. Regional maps could be produced based on 30 to 60 days of mission data, 
allowing investigation of temporal changes with monthly variation. 

This report lays out the methodology for this mapping analysis, providing the necessary 
elements of potential theory, orbital mechanics, and estimation theory. It is apparent that with 
0.5 p/sec SST data, centimeter accuracy geoid heights, with 50 to 100 km spatial resolution and 
monthly temporal resolution, can be recovered. Though, there are still details to work out, this 
objective can be met. 

The concept of a Satel 1 i te-to-Satel 1 ite Tracking (SST) mission is simply to calculate the 
gravitational force acting on a spacecraft from changes in its measured velocity. The satellite 
itself is the sensor, and its velocity, the observable. This is pictured notionally as follows. 
Consider a satellite in orbit around the earth approaching a region of excess mass. As the 
satellite approaches, it is accelerated toward the mass, and after passing the mass, it is 
decelerated. By measuring the time history of the velocity variation, an estimate of the 
magnitude and position of the mass excess can be deduced. Or course, the actual situation is 
much more complex for several reasons: the structure of the earth's mass distribution is very 
complicated, other forces act on the satellite, only one component of the satellite velocity is 
measured from another satellite, and the observations contain errors. In the SST concept, the 
second satellite could be very high, say in geosynchronous orbit or on another planet - the high- 
low configuration. Alternatively, the second satellite could be in the same low orbit, trailing the 
first low satellite - the low-low configuration. In the low-low case, the second satellite would 
experience similar velocity changes, but at a later time. 

SST tracking has been realized in a number of missions. The earliest example of 
Doppler tracking a satellite from a point not on the body being studied is the mapping of the 
lunar gravity field by means of lunar orbiters tracked from the earth [Muller and Sjogren , ( 1 968) 
]. Since this remarkable success, there have been notable Doppler tracking experiments for 
earth satellites: the tracking of Geos-3 from the ATS geosynchronous satellite. These are 
examples of the high-low configuration. There have since been numerous examples of tracking 
around Mars, Venus, and Eros: also examples of the high-low configuration. 

A number of analytical strategies have been employed. The most common is to treat the 
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SST data in the same way as other tracking data, and determine spherical harmonic coefficients 
of the body under consideration by variants of the method of least squares: an estimation 
problem in a finite-dimensional vector space of parameters. There are significant computational 
issues with global solutions [ Bettadpur et al., (1992) ]. In any case, these solutions would 
benefit from the high accuracy of the SST tracking data, but do not make any use of the unique 
geometry of the SST data. These ideas have been applied to estimating gravity anomalies 
[ Vonbun el al.J 1 978) ]. Alternatively, the analysis can be done as a mapping of the acceleration 
field at the satellite's altitude [Muller and Sjogren, ( 1 968) ]. Studies have also been done on 
mapping the acceleration of potential at the earth's surface [Colombo, (1981) ; Rummel, (1975) ; 
Rummel el al.,( 1 978) ]. Studies of recovery of fourier coefficients at satellite altitude have also 
been made [Kaula, (1982) ]. 

The analysis described in the sequel employs a generalized inverse theory (GIT) [Backus 
and Gilbert, ( 1 967) ; Backus and Gilbert, ( 1 968) ; Backus and Gilbert, ( 1 970) ; Backus , ( 1 968) ] 
to convert the observed velocity measurements to a mapping of surface values of geoid height. 
The ill-posed boundary- value problem and the unstable downward-continuation problem 
[Bullard and Cooper, ( 1 948a) ] are addressed in a way that the additional assumptions used to 
obtain the solutions are clearly identified and an error estimate is found. On the way to this end, 
a semi-analytical solution of the problem of calculating satellite perturbations caused by the 
anomalous potential will be formulated. Such a solution provides insight into the proposed 
measurement mapping, and significant computational economies. It gives a direct way to 
calculate analytically partial derivatives for the observable (velocity) with respect to the desired 
end product (the geoid). The analytical solution also allows the kernel function in the 
generalized inversion to be calculated directly. 
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II Background 


It is seen to be obligatory in this section to provide a prolix, far ranging, penetrating, and 
relevant review of the scientific objectives and uses of the analysis to be described in the sequel 
as well as a similar critical - in the classical sense - review of the past work. Since there are a 
number of carefully written reports available, [aww., (1979); anon., (1997); Kaula, (1970); 
Nerem et al., (1995 Aug 10)] this discussion will be limited to a few remarks on the method. 

Before, October 4, 1957, individual surface gravity measurements were made and 
essentially local gravity or geoid maps produced. Attempts to produce global gravity, or geoid, 
models from these data produced rather simple ones. Tracking of artificial earth satellites very 
soon produced significant improvement in knowledge of the global gravity field. The difference 
was that satellites provided a global sensor, that allowed measurement of large scale features of 
the geoid. In the succeeding years, analysis of tracking data has produced remarkable results, all 
based on using the integrating property of satellite motion. Combining the large scale 
information, from satellite tracking, with the small scale information from surface 
measurements -- surface gravity or more recently satellite altimetry — has produced combined 
models of remarkable resolution and accuracy. This has required development and refinement of 
algorithms and orbit computation methods for determination of spherical harmonic 
representations, which will be referred to as classical satellite geodesy.. 

Satellites provided two significant advantages. First, orbital motion samples large scale, 
long wavelength, features of the geoid. Second, satellites can quickly sample the geoid 
everywhere. The idea that a satellite samples the geoid everywhere was first dramatically 
shown, not for the earth, but for the moon by Muller and Sjorgen [ Muller and Sjogren, (1968) ] 
with observation of MASCON’s. Though the geometry is quite different - the lunar orbiter was 
observed from the earth where this analysis has two co-orbiting earth satellites, one tracking the 
other - the principle is the same: viz that instantaneous relative velocity changes are due to the 
difference in potential of the two satellites. So, to proceed in a rigorous way, one needs a 
mathematical formalism for: a) the potential mapping, b) the relation between potential and 
satellite velocity, and c) an estimation procedure. The estimation procedure chosen. 

Geophysical Inverse Theory Spectral Expansion Method (GITSEM), seems to make fewer 
assumptions about the form of the model, defining representation basis functions derived only 
from the observation geometry. It also provides a rigorous method of solving the downward 
continuation problem. The end result is intended to be a direct mapping of the satellite-to- 
satellite (SST) tracking measurement to geoid maps - the choice of geoid maps or gravity 
anomalies being quite arbitrary. The method can be used for any size map, including a global 
one, though in this case there seems to be little advantage of the method compared with a 
spherical harmonic expansion. In addition, a rigorous error estimate is also available. 

The view represented in this analysis, is to build on the results already available from 
classical satellite geodesy geopotential solutions. At present our knowledge of the low degree 
and order terms in the spherical harmonic model of the geopotential is quite good. We will 
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assume that up to a certain degree and order, a reference field can be adopted. In fact, the theory 
and implementation allows one to calculate local corrections to this model if needed. For 
example, in this analysis a reference field to degree and order 40 and another to degree and order 
90 were used. The derived map is then a correction to this reference model. 
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Ill Elements of Potential Theory 

The underlying mathematical framework of Geophysical Inverse Theory is based on 
potential theory. There are two facets of this discussion: 1 ) The solution of the boundary value 
problem in Geodesy, and 2) The downward continuation problem. There is an extensive 
literature on both of these topics {[Heiskanen and Moritz , (1967); Jekeli, (1981 ); Moritz, 
(1980,1989)]), and so we limit ourselves to a few remarks to define the notation and relations 
used in this analysis. 

From classical potential theory [ Kellog , (1953) ], the geopotential, T(r,<p,/,), - 
depending on the spherical coordinates r, the distance from the origin i.e. the center of mass of 
the earth, <p, the geocentric latitude measured from the equator, and k, the longitude measured 
from an arbitrary point on the equator [I A U, (1983) ] - is harmonic in free space, i.e. is a 
solution of Laplaces equation in three dimensional space exterior to all attracting masses, has 
continuous second derivatives, and is regular at r=°°. Therefore, T is analytic. There are a 
number of issues of principle concerning other masses (the earth’s atmosphere. The Sun, Moon, 
planets, galaxies etc.) that we treat with various approximations, viz, the atmosphere is modeled 
as an equivalent surface layer, the sun, moon, and planets are moved to °° by removing their tidal 
potential, and the remainder of the universe is regarded has having a negligible effect 
Therefore, given the values of the potential everywhere on a known surface enclosing all the 
earth’s attracting masses (the Dirichlet boundary value problem), T is uniquely determined 
everywhere outside that surface. 

If we chose the surface to be a sphere of radius R, centered at the earths center of mass, 
then T is given by Poissons Integral Formula [Heiskanen and Moritz, (1967), p35; Jeffreys and 
Jeffreys, ( 1 956,p22 1 )] 




/?(r 2 -/? 2 ) 
4/r 


/l 

J I 


T( R, (p , X ) cos (p d(pdX 


(3.1) 


/. =o <?=- n f r 

where, p, the distance from the integration point on the sphere (R,cp’,A.’) to the sample point 
(r,q>,A.) is given by 

p" = r 2 + R 2 - 2rR cos t// (3.2) 


and 

cos if/ = sin (p&mtp +cos ^cos (p cos( A — X ) i 


'N.B. We have written these expressions in terms of the latitude, (p, whereas they are 
often written in terms of the colatitide, 0=n/ 2-<p. In this report we will freely use both 
conventions, strictly adhering to this notation. 
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For example, R can be chosen to be: a) just below a satellite orbit, or b) small enough to just 
enclose all the attracting mass (6384403m). 


Poissons Integral Formula makes no assumption about the mathematical form of T, only 
that it be harmonic and regular at infinity. Flowever for both theoretical and practical 
applications, much use is made of spherical harmonics: solutions of Laplaces equation by 
separation of variables. Now Laplaces equation in spherical coordinates is: 

„ dV d 1 ^ d¥ 1 d 2 ^ 

-- tan^- 


2r- 


+ ■ 


= 0 


dr~ dr d(p~ d(p cos‘ (p dX 

Assume than one can write *F as: 

T=«(r)P(«7)r(A) 

Introducing the constants n(n+l ) and m 2 (m and n turn out to be integers) we find: 

r n 

*-(«-!) 


R(r) = 


P(^)=^(sin <p) 

the ubiquitous Legendre functions, and 
COSW/l 

sin mX 

Conventional usage is to call n the degree of the Legendre function and m the order of the 

cos mX\ 

Legendre function. The combination P nm (s\n (p ) 


r(A) = 


(3.4) 

(3.5) 

(3.6) 

(3.7) 

(3.8) 


sin mX 


are sometimes referred to as surface 


spherical harmonics and 


1 / 


. n* 1 




cos mX 
sin mX 


are sometimes referred to as solid 


spherical harmonics. It has also common to denote the surface harmonics as: 

cos mX 
sin mX 

The functions Y nm are orthogonal. That is when integrated over the sphere: 


(3.9) 


/^ lm (sin (p) cos mAP sr {sw (p)cosmXda = 0 

JJ P nm {s\n(p)smmXP sr (^(p)^f^Xda = 0 

P„ m (sin^>)cos/wA/ > tr (sin^))sin/w/k/cr = 0 

with 

da - cos <pd(pd A 


s * n,or,r * m 
s * n,or,r * m 
always 


(3.10) 


(3.11) 


It is convenient to define fully normalized Legendre functions /^ m (sin (p) such that the integral 
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over the sphere is: 



Osinp) 


cos mA' 
sin mA j 


d<7 - 4 n . 


(3.12) 


We have: 


cos mA 
sin mA 


= J( 2 ~ <^o)(2w + 1) ; ^^i /Usin <p) 1 


^„,(sinp) 

with the familiar Dirac Delta function <$, = 


(« + w)! 
0j« * m 
\\n-m 


cosmA 
sin w/l 


(3.13) 

Formulae for calculation of P nm (sin(p) 


can be found books on Mathematical Physics and Physical Geodesy. Great care must be given to 
accurate calculation for high degree and order functions and their derivatives. 

It is easy to see that outside a sphere of radius R, enclosing all attracting masses, the most 
general solution, that is regular at °°, can be written: 


GM 


rY— — - 

rll t /> -< sin «’>[ c nm cos mA f ^ nm sm R ^ * 

f n ' f ' 1 


»=0|fO 


(3.14) 


where G is Newtons Gravitational Constant (6.67423±0. 00009 x lO^cmVfgram sec 2 ), and M is 
the mass of the solid earth, ocean, and atmosphere, GM=3. 98600.441 77±0. 00000000 10 x 10 20 
cm 3 /sec 2 . Now if, as supposed with Poissons Integral formula, ¥ is known everywhere on the 
sphere r=R enclosing all attracting mass, then ¥(R,<p,X) can be represented in terms of a surface 

harmonic series given by (3.14). The determined geopotential coefficients — C nm , S nm - then 
provide the representation of ¥ everywhere outside R. 


The downward continuation error (DCE) is usually framed in terms of this 
representation. The DCE arises from the need to know ¥ at or near the earth’s surface, a point 
where r^R, and (3.14) diverges. In principle, one can obtain a series, using spherical surface 
harmonics, that represents ¥, provided one assumes knowledge of the attracting mass 
distribution above the sphere of radius r. 

GM * n — — 

4 '~(r,p,A) = Yj Z ^ sin r ) cos + sin ^l 0 - r - R (315) 

r n=0 m=0 


It must be emphasized that the resulting representation, ¥ has potential coefficients that 
depend on r- the distance from the earth center of mass to the desired point: e g. the earth 


surface or geoid - i.e. 


c_(r) 


r £ 


R [Jekeli, (1981) ,p36]. This is not a representation in solid 


spherical harmonics, and that the coefficients in (3. 14) and (3. 15) are not comparable for r<R 
Other than a representation of the potential for r<R, given knowledge of the coefficients 

it is unclear how to use this representation. 


In the practical application of (3. 14) to satellite orbits, the infinite sum is truncated to 
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terms n <n for two reasons, one practical - it is impossible to numerically perform an infinite 
sum - and the other fundamental. The measurement(s) to be analyzed are of finite accuracy, 
and there is a point where higher degree and order terms have no sensible affect, and are 
therefore ignored 2 * . Therefore, no matter what data analysis techniques are used, e g. satellite 
perturbation analysis or direct mapping of acceleration of potential difference, the resulting 
geopotential is smoothed, or averaged, usually characterized by the limiting degree of spherical 


00 


harmonic representation, i.e. y V n . 

T"(r,^,/l) = — ^ sin ?)p- cosm/l + S- sin w 4 R ~ r ~ 

r n= 0r-0^^ 

There is also the partial sum of (3.15), i.e.: 

V^ir, ip, A) = — X X ^m (sin ) COS mX + $m sin ^ r - R 

r n=0 m=0 


(3.16) 


(3.17) 


Now for r>R, a spherical harmonic coefficient representation of T has two sorts of 
errors: 1 ) errors of commission, i.e. errors in the coefficients themselves, and 2) errors of 
omission, i.e. errors due to ignoring the higher degree and order coefficients. If we let the true 

A A — 

geopotential be 4* , and the true geopotential truncated be 'F" we can define the error of 
commission as: 



_ _ ip" 


(3.18) 


and the error of omission as: 


£,(T) = T - 'F = ^ X I !jj Osin pfC cos mX + S ~ sin "4 R - r £ K 


•r* r» 


n- n* 1 rtr=0 


(3.19) 


Now the error of commission is compounded by the present practice of evaluating (3. 16), for 
r<R to obtain smoothed values at the earth’s surface [ Bullard and Cooper , ( 1 948b) ]. which is 
convergent everywhere. However, the errors in the coefficients are amplified by (R/r) n for each 
order. If we contemplate fields of degree 360, then the error amplification of the degree 360 
terms for a geopotential representation derived from data on a satellite at 450km altitude would 
be ((6378+450)/6378) 360 =4.6x 10 i ° For degree 180 the amplification is only 2. 14x1 0 4 . Errors of 
commission will grow very dramatically. Though one expects the higher degree and order 
coefficients to be smaller than the lower degree ones, (Kaula’s rule of thumb [Kaula, (1966) ]), 
they are also more difficult of determine, and consequently the uncertainties will be larger. 

Second, we have the representation error: the difference between (3.16) and (3. 1 7): 

(3.20) 


2 This assumes that the size of higher degree and order terms - as measured, for example, 

by the degree variances - continue to decrease, as models suggest [Kaula, (1966) ] 
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which depends on the degree and order of series truncation. In a careful and detailed analysis 
Jekeli, [Jekeli, (1981) ] has shown two critical results. Jekeli combines the error of 
representation (3.20) and the error of omission (3. 19 suitably generalized for r<R as 3.17) as the 
downward continuation error - the error of commission is not included. Jekeli, with n = 300 , 
finds (tables 3 and 4) that the DCE largest near the poles a(e t ) <0.090 mgals (0.290 mgals max) 
gravity anomaly and a(e r ) <0.042 cm (0. 14 cm max) geoid height. Second, by seeking 
anomalies averaged over a spherical cap of about 1 .4 degrees (tables 6 and 7), the DCE - again 
largest at the poles, is estimated to be a(e r ) <0 004 mgals (0.014 mgals max) gravity anomaly 
and a(e r ) <0.0020 cm (0.0066 cm max) geoid height. Therefore the conclusion of Jekeli 
[Jekeli, ( 1 98 1 ) , pi 27] “The downward continuation errors depicted in tables 3 through 7 are 
completely insignificant with respect to anticipated measurement accuracies of 1 mgal and 10 
cm in the gravity anomaly and geoid undulation, respectively.” And, “ .. the estimation of point 
or mean gravity anomalies and geoid undulations (height anomalies) using the outer series 
expansion to degree 300 anywhere on the earth’s surface is practically unaffected by the 
divergence of the total series.” This issue was discussed again [Wang, (1997) ] who confirms 
Jekeli’s analysis: viz “ .. , the method of smoothed analytical downward continuation can be 
used to determine the earth’s gravitational potential to any required accuracy.” 

The basic potential theory formalism has been presented, and the relevant formulae 
defined for use in this analysis. The downward continuation error of commission is controlled 
by suitable averaging, which is inevitable given a finite spacing of the data and a finite number 
of observations. Therefore we adopt the earth equatorial radius (R=6378137 m) for the 
reference sphere. We see that the essential problem of downward continuation is controlling the 
growth of numerical error of commission. We turn to data analysis methods that provide tools 
for this. 
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IV Orbital Theory 


In discussing the recovery of gravity anomalies, Ag, or geoid heights, N, with SST data, it 
is convenient to have analytical formulae relating satellite position and velocity (the state vector) 
to the desired quantity. Such formulae also enable the sensitivity, or partial derivatives of 
gravity anomalies and geoid heights to be more easily obtained. Otherwise, one must resort to 
costly numerical methods, with consequent loss of generality and insight. With this motivation, 
we proceed on two levels. First, some simple illustrative relations are obtained. Then a more 
complete treatment will be developed, for use with the actual data analysis. 

Consider, first, the theorem of conservation of energy for the satellite orbit. Of course, 
drag and radiation pressure perturbations would have to be taken into account if the following 
relationships were used for analysis of actual data. Let the total potential be represented by T, 
which for convenience, can be separated into a reference potential, p/r+U, and an anomalous 
potential, T, i.e. 

y = E + U + T ( 4 . 1 ) 

r 


In the conventional physical geodesy notation, U is the normal potential corresponding to that of 
a reference ellipsoid. However, here we prefer to view it as a reference, or adopted potential, 
with T being the remaining (unmodeled) part that we seek 3 . If we write the kinetic energy as 
'Arrw 2 , with the vector components of velocity along track (v u ), cross track (v w ), and radial (v s ), 
then 


2 



+ v. 2 +vJ) + H' 


= constant 


(4.2) 


For a satellite with small eccentricity, we can treat the along track velocity as the unperturbed 
velocity, v„, which gives v u =v 0 +5v u , v s =8 v s , and v w =8 v w . Therefore, to first order in small 
quantities, we have 


’Some confusion is bound to occur because the anomalous potential, or disturbing 
function, is generally denoted by R in celestial mechanics, and T in geodesy. Also the sign 
convention for potential in physics is reversed from that in celestial mechanics and geodesy. 


Here the force F = VHK = m — — . Finally, for convenience in this section, we refer to the 

dr 


product GM=p. 
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(4.3) 



This formula was first derived by Wolf [Wolf (1969) ] and subsequently used by many others. 
For small eccentricity satellites, equation (4.3) gives the change in along track velocity as a 
function of the potential in space (i.e. position and time) to about 10%, the errors arising from 
the change in radial distance owing to T interacting with g/r+U. Equation (4.3) is a direct 
mapping to the potential into the along track velocity, and could be used as a first approximation 
for inverting observations of 8v T to determine T. 

For a more complete theoiy, following the treatment in Brouwer and Clemence [Brouwer 
and Clemence , (1961)], using the potential, we can write the equations of motion as 


d 2 x x %U + T) 

~TT + ~ T 

dr r dx 

d 2 v v d(U + T) 

dr f 4 - 

drz x d(U+T) 

dr r c& 


(4.4) 


where p=GM=3.986xl0 20 cm 3 sec' 2 and 

r 2 = x 2 + y 2 + z 2 (4.5) 


If the coordinates are r, u, and <p, where the longitude from the equator crossing in the orbit plane 
is u=f+o), (true anomaly, f, plus argument of perigee, w), and <p is the latitude 4 then the 
equivalent differential equations of motion are: 


d'r 

~-r cos' cp\ 


dt 


f-V-r 

dtp 

[dt) 

dt 


H %U + T) 


+ — - 




d ( 2 du\ c\U + T) 

— r cos' (p — = 

V dt) ck 


dt 


d ( dqr 
— r' - J1 
dt V dt J 


+ r sin (p cos (p 


du) d(JJ + T) 


KdtJ 


dcp 


(4.6) 


4 The derivation is general. Below, the orbit plane will be used as a reference, and 
perturbations in along track (5u), cross track (8w) and radius (8r) will be developed. In this case 
the latitude, cp, will become the cross track, w, component. 
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To obtain an expression for the potential, U+T, as a function of position (and time), we 
begin by assuming that the potential on the earth’s surface at r=R is (U+T) lerr (R,<j>A), where R is 
the radius of the earth and k is the longitude To express T in an inertial system appropriate for 

equation (4.4) or (4.6), we have ( U + T) lner1 {R,<p,A-0(t)), where #(7) is the sidereal angle at the 

time in question. For the moment we ignore the effects of a moving equator [Gaposchkm, 
(1973)]. 

Following the discussion in section III, since (U+TXR,<p,^) given on a sphere is the 
Dirichlet boundary condition for Laplace’s equation external to the sphere r=R, we can use the 
basic results from potential theory to obtain T at any point in space outside this sphere, see 
section III. Therefore we can expand T in terms of orthonormal base functions (associated 
Legendre polynomials) as 

n- x mn 

11 

n-2ffh0 

which, using the properties of solutions of Laplace’s equation in spherical harmonics, can be 
upward continued as 


p (sin fit) 
nm 


C cos mX + S sinw/l 
nm nm 


(4.7) 


(U + T)(R,<pJ) = ^- 


(U + T)(r,<p, A) = 


GM 


fi-x n ( n\ n 

XX* 

n= 2 m=0 v ' 


F C cos/wi + S sin mX 


p (sinc>H 

nm l nm nm 


(4.8) 


Alternatively, since U+T is harmonic in space, we can use Poission’s integral to obtain 


(U + T)ir ^ X)= BrlsJll 2 \ "j (u±lIM^l c "tw ( 4 . 9 ) 


4 7t 


\ \ 

/. '= 0 <p'~- t/2 


where p is the distance between the integration point R,<p\A.’ and the sample point r,<p,X; it is 
often written as 


p 2 - r 1 + R 2 - IrRcosy/ (4 10) 

vp being the central angle, which can be expressed 

cos^ = sm^cos^>'|-i-cos^cos^ , cos(/l-/l'). (4.11) 


Now, assume the motion of the satellite is given in two parts: i.e. 
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(4.12) 


x = x Q + dx 

y = y 0 + fy 

z = z 0 + dz 

where x^y,,, and z 0 satisfy 


d x a x 0 

-TT+P- f 

dt r a 

d 2 y y 

— + Li — 

j.l r- 3 


+ M~f = 


and r=r +5r. We also have 


, nj %y + T). ftJJ + T), <%U + T), 

d(U + T) = — -dx + — -dv + — -dz 

dx dy dz 

.... c\u + t), a u + t ), au+T), 

d(U + T) = — -dr + — -du + — -d<p 

fr du dq> 


au+T) au+r ) au+T) au+r) 

r— - = x — -+y— - + z — (4. 

a dx ' dx dx 

Now, multiplying (4.4) by 2dx,2dy, and 2dz respectively, adding, and integrating the 
result we have 


dzX 2 // , n 


+ ^- = ljd(U + T) 


where p/a is an arbitrary constant of integration. This is chosen such that 


jy J*m) -lK + tL = 2 tou 

v dt ) v dt y r a J 
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is satisfied. In terms of polar coordinates (4. 16) becomes 


dr^ 2 

dt ) 


+ r 2 cos 2 <p\ 


r duY ■>( d(p Y 2 n fi 


VdtJ 


+ r‘ 


\dt) 


+ 


a 


Again, multiplying (4.4) by x,y,and z respectively, adding 


(U + T). 

sing (4.15) we have 


d 2 x d 2 y d 2 z n d(U + T) 

x — r + ^-rf + z ~r- r + - = r-b — — - 

dr dr dr r cr 


or in polar coordinates 

d 2 r , , (du'Y 

r r~ cos" (p — 

dt \dtJ 


2 ( dtp V fj c\b 

- r — + — = r 

\ dt J r t 


Adding (4. 18) and (4.20) we have 

L^l-t + t = 2 \d(U + T) + r 
2 dr r a J 


c\U + T) 
dt 


Now using r=r 0 +8r, and subtract the reference orbit 


1 fu [ // 

2 dt 2 r . a 


= 2 JdU 


+ r 


dU_ 


to first order in small quantities 8r and T we have 


d 2 {r G Sr) fi r , T , <% . n 

-^- + ir>-2jrfr+r.-+a-— -ft 


dr 




We also have 


(4.18) 


(4.19) 


(4.20) 


r id 


(4.22) 


(4.23) 
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(4.24) 


d 2 5x 

+ £. 

dt 2 


d 2 Sy 

+ JL, 

dt 2 

V 

d 2 dz 

+ JL 

dt 1 

y 


( 


1 1 


( 


1 1 


fJ r. dr 

dx yr n 

dT_ 

dy ' Vr- 

dT f i 1 

~ At ^ ^ 3 3 

dz \r n r 


a 

Qy 

Q: 


Now, a solution to (4.23) and (4 . due to Hansen, can be found {[Brouwer and 

Clemence , (1961)] chapter XIII) for the mutations in radial (8r), along track (5u) and cross 
track (5w) and can be written: 


Sr=\Xsm[f-fW 

Su= j[Jiy/-2f]rf/ (4.25) 

dw= Jz sin[/ -/']#' 

where 


r 3 r 3 

X = — O r = — 

VP~~ VP 


cT „dU 
2 \dT + r — + Sr 
J dr 


cY 


Y = 

Z = 


r_dT 
fjp du 
d_ dT_ 

fjp dw 


(4.26) 


where p=GM=n 0 2 a 0 3 and p=a( 1 -e 2 ). Here, the integration variable is the true anomaly, f, but for 
practical numerical implementation, we will use time as the variable of integration, using 



(4.27) 


In equation (4.26) we interpret dT as 


dT ,dTj ^dT . 

dT = — dr 4 du H dw 

dr du dw 


(4.28) 
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i.e. the integral along the satellite path. Therefore fdT = A7\ the difference in potential from 
a reference point on the orbit to the point of interest. 

Now, we measure relative velocity, and can obtain the perturbation in velocity - in 
radial, across track, and along track direction - as: 


5Vr eh 5 *' (1 - e 2 ) m 


'"-Vjj X cos( f ~ f )df'-X sin(/ - /„ ) 


\r 


Jo 

n 


SVu r dt 3 " \\-e 2 ) m 


'a'* 


\rJ 


■\Ydf- 


2 dr 


fo 


(4.29) 


o d s n 

^ =r j t Sv= j^r 


^v' / - 


\r. 


)z cos(/ - /' )df'-Z sin(/ - /„ ) 


Jo 


For illustrative purposes, these expressions can be simplified for small eccentricity. Of course, 
in all numerical calculation, the full expressions (4.25) and (4.29) are used. The simplified 
expressions are: 


X = 


n 2 a 


rrr 

2 f dT + a — + Sr 

i dr 

n 2 a dM 

Z = -L— 

n'a ck 


dU 


(4.30) 


Now 


* ■ -h\ XMf - f)df '= sm</ ' r w 


(4.31) 


and 
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(4.32) 


Sv r = «{ jx cos( M - M' )dM'-X sin( M- MJ J 
<Sv u =/7{j}™-2<5r} 

. = ra { fz cos( M - M' )dM'-Z sin( M - M 0 )} 

Equation (4.3 1 ) defines 5r in terms of (5r/a)( <3U/dr), and can be solved by iteration, as (eU/cT) is 
of the order 10' 3 . Now we can rederive equation (4.3) using the second parts from equations 
(4.30) and (4.32) 

dv u = n\ a ■ - 1 - T- 25r \ - — T - 2 nSr (4.33) 

V n~a‘ J na 

where we can make the identification v 0 =na. Now the second term can be written 


* GM GM 

nor = t-ot = -or - 

na naa" 


GM_&__]_d_( Gh£ 

v„ a 2 v„ dr V r > 


dr 


(4.34) 


which is simply the change in velocity due to the perturbation, 6r, with the central force term, as 
pointed out earlier 


We can now formulate the observable, A p . We 
have the satellites, P and Q, with position r p ,r Q and 
v p . v a . The relative position between P and Q is 


P = r P ~r Q , 

the distance between them is 


(4.35) 


P=y[p 7 P, 

and the relative velocity is 


(4.36) 


P-P-(v r -v L> ). (437) 

Now, we suppose that we have a reference orbit, 

r Po , r Qo » v Po , v Qo , and the true orbit, r p , r Q , v p , v Q 



The difference, or perturbation, to be modeled is dr p , Sr^ , Sv p , SVq . Now consider the residual. 


or difference between the observed, assumed to be the true, relative velocity, and the relative 
velocity computed from the reference orbit 


Ap=p-p o =p-p a *( v Po -v Qo ) 


(4.38) 
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To define the algorithm for calculation of the perturbed position and velocity, F, v , from the 
unperturbed position and velocity, F 0 ,v o , and the along track, 5u, across track, 5w, and radial 
perturbations, 5r, (Figure 4. 1 ). 


We use the unit vectors 


f = r 


r -r 


w = 


r xv. 


JTr x v) (r x v) 




u = w x r 


(4.39) 


for the radial, cross track and along track directions. More precisely, the along track unit vector 
is normal to the radius vector, and not along the velocity vector, to be consistent with the 
definition of the along track perturbation. So we can write 


~du 

JK 

/v — 

S — 

dt 

= u r v ] 

— U * V 

o v c 

e dw 


/A — 

S — 

dt 

= w, - V, 


„dr 

A — 


5 — 

dt 

= r r v i 



These relations lead to the perturbed velocity 


f A A A 1 

v . = h w i r .j 


u 

o 


VV 


f. 






i] 


. du 

~dt 

dw 

dt 
, dr 

~dt 


(4.40) 


(4.41) 


Now, with two satellites (P,Q) at positions x P and Xq, the unperturbed satellite to satellite 
range p 0 is 

Po = yj{xro-x Q o)-{xpo -*Qo) ( 4 42 > 

from which we can obtain the unperturbed range rate 
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(4.42) 



The corresponding perturbed range rate is 

P\ = P\ •(*/>> - ) (4-42) 

If we define 

A x, = [«v, w AV r Xj ] (4.43) 

the perturbation in range rate, A/7 = p x — p o , is 

A/? = • { A Pl A p J Po + A P] [duwr p ] - A gx A' Qo jc Qo - A Ql [Suwr Q ]} - p Q • {jc Po - * Qo } 


With this formalism, we give an example of this calculation. The full Cartesian 
equations of motion for perturbed (true) and unperturbed (reference) orbits are integrated, in this 
case with a 4 th order Runge-Kutta integrator with variable step size error control, set at 10' 16 , 
using 20 digit precision calculation. Both models include the quadrapole moment, J2. The 
twenty four equations of motion are integrated. The quadrature of (4.25) and (4.29) is done 
following the presumptive data to be obtained from the GRACE mission. We assume equal time 
spacing of the data, nominally 10 seconds, and use a modified Simpsons rule [Hamming, (1973) 

] quadrature formula. Using the unperturbed position and velocity, the six perturbation 
equations (4.25) and (4.29) are integrated. In fact the time spacing for the Simpsons rule was set 
to 1 sec, 2 sec, 5 sec, and 10 sec. The 10 sec quadrature interval had some error build up: for 
actual data analysis, the perturbations should be calculated with a 5 sec or smaller time step. 

This is not difficult, as the reference orbit can be obtained at any time resolution. 

For illustrative purposes, two perturbing potential cases will be shown. The first, is a 
single geoid height anomaly in the subsatellite path, and the second is with two geoid height 
anomalies separated by 1 degree, =110 km, of equal size and opposite sign: a classical dipole. 
These will illustrate some important properties of the perturbations, to be used in developing the 
data analysis methods. Table 4. 1 gives the initial orbit elements of the trajectories (the cartesian 
coordinates are used), and table 4.2 defines the anomalous potential. 
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Table 4.1 Initial Orbit Elements 


Epoch (MJD) 


x (cm) 


dx/dt (cm/ sec) 


y (cm) 


dy/dt (cm/sec) 


z (cm) 


dz/dt (cm/ sec) 


a (cm) 


e 


I (deg) 


M (deg) 


O (deg) 


to (deg) 


Satellite P 

Satellite Q 

51478.000 

51478.000 

6.774373333781 1697E+08 

6.774373333781 1697E+08 

-2.2620699340959200E+04 

2.2620699340959200E+04 

1 .04627593 19408060E+06 

- 1 .04627593 1 9408060E+06 

4.0 121 0871 68058562E+04 

4.0121 087 1 68058562E+04 

1 .9964 1 34070442688E+07 

-1 .9964 1 34070442688E+07 

7.6555594831 3981 76E+05 

7.6555594831 398 176E+05 

6.778E+08 

6.778E+08 

0.0001 

0.0001 

87.0 

87.0 

1.69 

-1.69 

0.0 

0.0 

0.0 

0.0 


Table 4.2 Anomalous Potential 

Case 

<l> (deg) 

X- (deg) 

H (cm) 

Size (km) 

I 

45.0 

325.0 

400.0 

110.0 

II 

45.0 

325.0 

400.0 

110.0 

II 

46.0 

325.1 

-400.0 

110.0 
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perturbed position and velocity be x px , v Pl . 
So in figure 4.2 we show the change in 
velocity 5 , 

Sv P = v Pl - v Po for satellite P for case I, and 
figure 4.3 shows the change in radius vector, 

8r p = r„. - r Po for satellite P. These are 
obtainea from numerical integration of the 
equations of motion. If the analytical theory 

for 5v p ,8r p were plotted on the same graph, 
they would overlay the numerical result. 

Figure 4.4 shows the difference in the 

numerical and analytical theory for 8v p ,8r p . 
The rms difference between the numerical and 





theoretical perturbations for the two curves is 2.4xl0' 5 (cm/sec) = 0.24 (p/sec) and 0.019 cm. 
The long period structure in Figure 4 4 is an unmodeled interaction perturbation between 8r and 
J2. With J2=O.0, the numerical and analytical theory agree to 10' 10 . The high frequency 5v 
noise is due to rounding error in data input to the plotting program, and not an orbital effect. 

Figures 4. 5,4. 6, and 4.7 are the same quantities for case II, the dipole anomaly. The very 
small 8r perturbation (0. 1 8 cm maximum) results in negligible J2 interaction and long period 
terms. The reduced Sr also results in smaller contribution to 5v. For case II, the rms difference 


5 Strictly speaking we are plotting speed. Velocity is a vector quantity, and we measure 
and analyze a scalar quantity, the speed. However, common usage is to refer to this scalar 
quantity as velocity. In this discussion, vector quantities are always given as x , or x for a unit 
vector, and scalar quantities as x. 


- 23 - 


i 







between the numerical and theoretical 
perturbations for the two curves is 2.7x1 O' 6 
cm/sec and 0.0006 cm. N.B. that 0.0006 cm 
gives 6.4x1 0‘ 7 cm/sec contribution to the 
velocity. 

We now illustrate the difference in 
velocity of two satellites, see table 4. 1 , 
separated by 400 km, at an altitude of 400 km. 
In figure 4.8 we show the difference in velocity 


° Geoid Recovery 

Cm I 

ftiagi* 1 <bynv 44W cm Geoid Height Aiotnoh 


if\ ' 

I \ dT T(T>-T(Q> 


w*dr-n(tlr F -AQ) 



for case I. Note that here we plot the difference in absolute velocity, not the rate of change in 
distance between P and Q. See section V for further discussion of this point. We also plot the 
difference in anomalous potential, AT, between P and Q, and the contribution 
ofnAr- n(Sr p - 8r^) to the difference in velocity. Recall that Av * AT + nAr The root 

mean square, rms, difference between the numerical and theoretical perturbations are 2. 5x1 O' 4 
cm/sec in differential velocity and 0.0024 cm in differential radius (2.6x1 0 -6 cm/sec as 
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velocity). Figures 4.9 shows the same quantities for Case II. For Case II, the rms difference 
between the numerical and theoretical perturbations are 2.9xl0' 5 cm/sec in differential velocity 
and 0.0002 cm in differential radius (2.3x1 O' 7 cm/sec as velocity). 

Jekeli [Jekeli, (1999) ] has studied this issue using the energy conservation theorem with 
quite different objectives. He presents an expression for the corrections to Av to obtain AT. In 
this case he proposes to measure these corrections “in situ” with the GPS receivers on the 
GRACE satellites. However, the requirements for the “in situ” measurements exceed the 
present GPS capability. 

The simple simulation shown here demonstrates a number of facts. 

1 ) First we see the, well known, sensitivity of relative velocity to the anomalous 
geopotential. Here we have postulated a rather large lxl degree 400 cm geoid anomaly. 

2 ) The significant indirect effect of the radial perturbation on the relative velocity. To 
make best use of the relative velocity measurement, one must treat this indirect effect. The 
radial perturbation also has a significantly different time history (or fourier spectrum), with the 
presence of long period effects. This should compared with the potential difference, AT, 
contribution, which is more than 1 0 times larger, and is a local effect confined to the dimension 
of the anomaly. This will be discussed further in section V. 

3) The efficacy of these formulae in calculating the perturbations, given the anomalous 
potential. One use of these equations could be to correct the observed velocity for the effects of 
an anomalous potential in an iterative procedure. The relative efficiency of such a calculation, 
requiring quadrature over the total anomalous potential field, compared with use of a spherical 
harmonic representation would have to be investigated. 
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V Satellite to Satellite Tracking 


To define the algorithm for calculation of the perturbed position, x„ and velocity, v,, 
from the unperturbed position, x ( „ and velocity, v„, and the along track, du, across track, dw, and 
radial perturbations, dr, (Figure 1 ). 


If we use the unit vectors 


- r 


r r 


w - r x v. 


AAA 

u = w x r 


yfcr X v)-(r X v) 



for the radial, cross track and along track 
directions. More precisely, the along track unit 

vector is normal to the radius vector, and not along the velocity vector, to be consistent with the 
definition of the along track perturbation. 

More precisely, the along track unit vector is normal to the radius vector, and not along the 
velocity vector, to be consistent with the definition of the along track perturbation. So we can 
write 


.. du 

A 


5 — 
dt 

= «i • v. 


~dw 

A — 


8 — 
dt 

= w, • v, 

- w 

x dr 



8 — 
dt 

= • v, 

~K 


v 


O 


* V 


o 


V 


o 


These relations lead to the perturbed velocity 


v i=["i w, r,] 


U 


V. 

K 




w, 


n] 


, -.du 

o — 

dt 

c dw 

s Yt 

s* 

dt 


( 5 . 1 ) 


( 5 . 2 ) 
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Now, with two satellites (P,Q) at positions x P and Xq, the unperturbed satellite to satellite 
range p 0 is 


Po=tJ{xpo-Xqo)‘( x Po-X Q o ) 

from which we can obtain the unperturbed range rate 

a _ 

dt 


Po=Po'{ X Po- X Qo) 


The corresponding perturbed range rate is 

P = A •(*/». "%) 


If we define 


A \ •,=[“» r Xi ] 


the perturbation in range rate, Ap = p — p a , is 


(5.3) 


(5.4) 


(5.5) 


(5.6) 


Ap — p x -^Ap^ApoX po + Ap^duwrp] A () 1 Aq 0 Xq 0 j| P 0 ’\xp„ -^oo} 

(5.7) 


The basic observable desired is the difference in velocity, i.e. the intersatellite range rate 


is 


Av = y[x~p • X p - Jxq • X Q 


(5.7) 


where the SST measurement is give by (5.5). If we assume that the perturbation is along the 
velocity vector, then these are related by 

Ap = QOS(0 p )y[¥ p • Xp - COS(0 Q )yjx Q • Xq = COS(0 p )v Q - COS(0 Q )v c; (5.8) 

where 
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( 5 . 9 ) 



where 6 = (G p + 6 0 ) I 2 . This correction was applied in the following simulations, resulting 
in small improvements consistent with the small size of 0. 
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VI Geophysical Inverse Theory 


The two issues concerning the fundamental process of converting range-rate 
measurements at satellite altitudes to geoid heights at the earth’s surface are: a) the solution of a 
boundary value problem in potential theory with incomplete data on an undefined surface, and 
b) the downward continuation of the potential to the earth’s surface in a manner that keeps the 
errors within bounds. A related issue concerns the regularization of the data so as to account for 
the mass between the geoid and the external boundary'. 

Geophysical Inverse Theory (GIT) allows one to address these issues in one step.. The 
GIT motivation and theoretical framework was first introduced by [Backus and Gilbert, ( 1 967); 
Backus and Gilbert, ( 1 968); Backus and Gilbert , ( 1 970)] and has had wide application. A good 
source for further developments and survey of applications and results, as well as the approach 
used here, the Spectral Expansion Method (SEM) can be found in [Parker, { 1977) ] and [Parker, 
(1994) ]. Therefore, we provide a brief discussion of the approach, and some details of the 
implementation for the problem addressed in this study. 

One central idea in GIT is to seek a representation of some geophysical quantity with a 
continuous function. The determination of a model with a finite number of parameters is a 
problem in statistics, e g. least squares and is not Inverse Theory. In practice, for this continuous 
function, with a finite number of errorless observations, an infinite number of models can fit the 
data exactly. For selection of one of these models, GIT provides methods for finding the model 
that minimizes a norm, i.e. in some sense the model minimizes some property of the solution. 
This is also generalized to data with errors when one would not want to fit the data exactly 
anyway. There are many possible norms to chose: study of this is beyond the scope of our 
discussion. The norm, used here is to find the solution that minimizes the square of model 
integrated over the sphere, i.e. the so called L2 minimization. We will also introduce a 
seminorm minimization, introducing arbitrary functions into the solution space that will allow 
specific model parameters to be determined. Norm minimization is also a basic concept in Least 
Squares Collocation [Moritz, (1978) ], and some comparison of the two methods is given in 
[Parker, (1994) ]. 

The Geophysical Inverse Theory Spectral Expansion Method provides a means to 
incorporate error analysis, an arbitrary amount of smoothing, and a direct mapping of the 
observable - in this instance range rate, 5v, converted to potential difference - to geoid height, 
N. As discussed in Section III, we can adopt a reference potential, U, including the principal 
oblateness and other low degree and order reference potential terms that we can assume known 
with sufficient accuracy. In addition we can safely assume a sphere of radius a e =6378. 137 km as 
the boundary surface for the boundary value problem. Finally we can assume the effects of mass 
between this reference surface and the observation are negligible. Furthermore, as discussed in 
section IV, the measurement of 8v between satellites P and Q, can be processed in a precision 
orbit determination calculation using the reference potential, U, to obtain residuals, Av. These 
residuals represent three factors, a) the potential difference, AT, between the satellites P and Q, 
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b) a perturbation, n8r, due to AT interacting with GM/r, and c) other orbital errors. By 
correcting for n5r, and filtering the other orbital errors, the residual, Av, gives an estimate of 
AT=T(P)-T(Q). 

In this work we consider the fundamental desiderata to be the geoid height, N, often 
referred to as geoid undulation, from the reference surface. To relate geoid height to 
geopotential we use Bruns formula [Heiskanen and Moritz, ( 1 967), page 85)] 



where y is gravity on the reference ellipsoid, strictly in this application it is gravity on the 
reference surface. Since we are concerned with small corrections to the geopotential, we use 
7=978.0327 (cm/sec 2 ) 6 to scale the geoid model, N, to geopotential, T. 

The Poisson Integral Formula (3.1) and Bruns Formula (6. 1 ) allows us to solve the 
forward problem, i.e. given a geoid model on a reference sphere, compute the potential at a point 
in space, P,. outside the reference sphere This is of the form: 

T(P,)= \g(P,S)T(S)da s =y\g{P l ,S)N(S)d<7, (6.2) 

cr a 

The P, is the position of the observation, S is on the reference sphere and do s is the surface 
element of integration on that sphere. The estimate of AT=T(P)-T(Q) is then written 

mP„Q t ) = r f[g(P„S)~ gtQ„S)]mS)d<r s =r \G(P„Q l ,S)N(S)d<r s (6.3) 

a o 

Anticipating the treatment of observation errors we can by use: 

A T /° 

and G, - (6.4) 

cr, cr, 

where o, is the standard error of the observation AT , . 

Now, the functional (6.3) is linear in the unknown, N, and therefore we have a linear 
inverse problem. Assuming P, does not lie on the reference sphere, (6.3) is also finite. With n 
values, or observations, we have AT i = 1,2, ... n there are n kemals in (6.3) 

i 

G,(S) = G(P„Q„S) / = U...« («•» 

Which are n functions of position, S. With a little thought one can see that these are linearly 

6 Geodetic convention defines the unit of 1 .0 gal as 1 .0 cm/sec 2 . 
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independent, i.e. G, is linearly independent from Gj if (P„Qj) is not the same as (Pj,Qj) 
Therefore they could form a set of basis functions to represent the geoid height on the sphere, 
i.e. 


N(S) = Y i a,G,(S) (6.6) 

7=1 

One can show that this has the minimum norm [Parker, ( 1994),Chapter 1], Substitution of (6.6) 
in (6.3) resulting in the matrix equation for a, . 

/[ r »] bH A? :] (67) 

where T is a symmetric matrix, called the Gram matix, with elements 

r „ = r ,. = jG,(S)G J (s)da !S (6.8) 

a 

This system of equations can be solved for dj . 

Now, the spectral expansion method, SEM, is introduced at this point, leading to great 
computational simplification and functional utility. The positive definite matrix T can be 
diagonalized with an orthogonal, eigenvector, matrix 0 with 

©T© = A (6.9) 

where 

A = diag{ Z l ,A 7 ,Z 3 ,....,A n ),A l > X, > A? > X n >0. (6.10) 

The matrix of eigenvectors, 0, has the usual properties 

©'© = ©©'=/ ( 6 . 11 ) 

the unit matrix I=diag(l,l,l,....l). The eigenvalues, ^ ,are often referred to as the spectrum of 
the problem. This diagonal decomposition can be accomplished by methods such as described 
in [Golub and Reinsch, (1971 ) ]. We can now define another set of orthogonal base functions 7 


7 The association of l/X^ with ^ follows [ Parker , (1994) ]. This change from associating 
1/Aj with both a, and , as in [Parker, (1977) ] is quite arbitrary, and is used to facilitate the 
seminorm minimization outlined below. 
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( 6 - 12 ) 


V',(/ > )=Z©;,G y (/ > ) 

J 

and write for the computed potential 

r</>) = Za,y,(P) 

I 

where 

a = 0 a 

1 Lmt 7 ' J 


Substituting (6.9) and (6. 14) in (6.7) we find 

a , = — V© /l -A7’ / 


(6.13) 


(6.14) 


(6.15) 


The a, are uncorrelated, i.e. are statistically independent, with the standard error of each 
coefficient 1/X, , i.e. 

or, =1/4 (6.16) 

These basis functions, (6. 13) become more oscillatory as i increases, i.e. as X, decreases, 
and approaches +0, see below for some examples. Consequently, the standard error of the 
coefficients increases as i increases. The power of the SEM lies in using Xj to select the desired 
solution. For example, by eliminating basis functions with high frequency oscillations - i.e. 
those functions with eigenvalues less than some minimum eigenvalue, - one can obtain a 
smoothed solution. Of course, not using the complete eigenfunction expansion, the model, 

(6. 13) would not fit the data exactly. An alternative would be to add a constant, C, to each Xj in 
(6. 15), which has the effect of significantly reducing the effects of small eigenvalues on the 
model. So, in matrix notation we have 

[fl]=[A+C/[ 1 [e'|AT] (6 17) 

In the case were there are observation errors, using (6.4) in (6.7) and (6.8) leads to a new 
solution if the observation uncertainties are different for each observation. Otherwise, the 
standard deviation of each observation cancels, and the result is the same. This leads to 
considerable computational simplification in this analysis, as we can compute (6.7) once for a 
given data set, and analyze the effects of random errors by scaling T, and adding noise to AT. 
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In the case were there are observation errors we can compute a statistic for the data set 


Z 2 =’£(*■?; ?f=n (6.18) 

1=1 

This can be used as the criterion for selecting or C. Having selected ‘k mm or C one now has 
the formal error of commission in the model 

<6I9) 

/=! A i + L 

where the obvious functions are removed based on selection of . 

The orthogonal functions (6. 12) used to represent the geoid are finite in number, and 
therefore define a limited functional representation of the geoid. For example they do not 
provide the same representation as the same number of spherical harmonics. The spirit of GIT is 
to let the data distribution define the representation. As an aside, it is the distribution of the data 
point positions (P,Q) that define the Gram matrix, (6.7), and not the actual data (AT), hence the 
representation. 

There may be functions that, a priori, we wish to include in the representation. For 
example, a constant, low degree polynomial, or spherical harmonic, that would not be included 
in the basis (6.12). For example, if the reference potential field, U, has a long wavelength error 
that could be locally modeled by a low degree bivariate polynomial. By including such 
functions, one cannot obtain a minimum norm, but can achieve a subnorm minimization. We 
can accomplish this expansion of the basis functions as follows. Let us include in the 
representation, a linear function: 


nA 

h(S) = Y.WS) ( 6 . 20 ) 

Jt= 1 

where the b k are to be determined from the data. The model becomes 


n nK 

V(S) = £>,«/, OS) +2>A(S) (6.21) 

/= i *=] 

Assuming that h k and y, are independent function of position, S, making the same substitutions 
as before, we have the system of equations to solve 
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( 6 . 22 ) 


"a + C7 

A 


a 


0'A T 

A‘ 

0 


b_ 


0 


where the elements of the rectangular (n x nK) matrix are 


4 = jA,(5)G/S)rf CT 

£7 

and 

[A]=[e'\A] 


(6.23) 


(6.24) 


As with (6.4) each element of (6.24) should be divided by the observation uncertainty when 
treating data containing errors. This system of equations has the solution 



and 


a'(a+ci)~'a 


^'[a+c/] '©'[at 


[fl]=[A+c/]‘ l [e'[Ar]-i<i 


(6.25) 


(6.26) 


In the application here, we intend to use low degree polynomials, centered in the region of 
interest (q> 0 ,A, 0 ). These polynomial coefficients are used to account for errors in the reference 
potential model, U. The errors are assumed to be long wavelength variations, and can be 
suitably modeled with low degree polynomials in longitude and latitude. The polynomial 
variables would be <p-<p 0 and X-X Q . Since the region of interest is of the order of a few degrees, 
(<0.2 radians) we choose the variables 


£=sin {(p-(p 0 ) 
rj - sin(A-A 0 ) 


(6.27) 


with the 21 polynomials 1 ,^,rj,^ 2 ,^rj,r| 2 ,...,£ > ri 4 ,ri 5 . This removes questions of definition in 
computing (6.23). 
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VII 


Simulations 


VII. 1 Introduction 

In previous sections, necessary elements of. Potential Theory (III), Orbit Theory (IV), and 
Geophysical Inverse Theory (VI) to analyze satellite to satellite tracking data for geoid 
determination have been reviewed. Now, some numerical experiments will be described that 
illustrate how these elements can be combined, and what sort of results may be expected. These 
numerical experiments will be progressively more complex and complete, in four stages. First 
will explore the nature of the eigenvalues and eigenfunctions that come from the singular value 
decomposition. This requires combining the elements of Potential Theory with Geophysical 
Inverse Theory. As seen in section VI, these results depend only on the physical geometry of the 
measurement, and are independent of 
the measurement itself. Second, some 
properties of the solution assuming a 
direct measurement of potential 
difference: accuracy, resolution, error 
propagation, and sensitivity can be 
shown. For these analyses, idealized 
geoid anomalies will be used: blocks 
with a dimension of 1 .0 degree and 5.0 
degrees will be combined. Third, an 
orbital simulation using a small number 
of idealized blocks is done, to illustrate 
geoid recovery using SST 
measurements. Finally a number of 
orbital simulations will be offered, with 
increasing complexity in the desired 
geoid, and other orbital effects. For 
these the geoid model will be the 
EGM96 [Lemome and al., (1998) ] 
geopotential model. EGM96 is a 
complete model to degree and order 
360. For the illustrative proposes of these demonstrations, the EGM96 model will be truncated, 
at degree 70, 90 and 1 80. In addition other physical forces will be included: lunar and solar 
forces, ocean and solid body tides, atmospheric drag, solar radiation pressure, and earth albedo 
pressure. The focus of these last simulations will be to demonstrate the recovery of geoid 
information from satellite to satellite range rate, which involves correcting for the indirect 
effects (equation 4.33) in the presence of other force model errors, that can only approximately 
accounted for in the orbit fit. Exploration of the degradation of the recovered geoid in the 
presence of these other forces, is a much larger effort than was addressed here, but some limited 
cases will be shown. 
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VII. 2 Geometry of GIT 

A set of data is generated with the parameters given in Table VTI_2_ 1 . 


Table VII 21 Parameters for Simulation 

Area Size 

20°x20° 

Geoid anomaly Size 

5°x5° 

Geoid anomaly amplitude 

8 cm 

Satellite Altitude 

300 km - 450 km 

Satellite Separation 

1°- 3° (110 km -330 km) 

Number of Data Points 

1000 

Distribution of positions 

Random 

Distribution of orientation 

Random 


Figure VIII is a cartoon showing this configuration of data. Figures VII_2_a,b,c, and d shows 
the eigenvalues for this system, often referred to as the spectrum of the system, for various 
combinations of orbit configuration. The enormous range, 10 14 , of the eigenvalues is evident. 
Note that the determinant of the system, the product of the eigenvalues, cannot be computed 
directly, as it is smaller than the minimum number represented in the computer (3.4x1 O' 1932 ). Of 
course the logarithm can be computed 
(-7075.2963). The Singular Value 
Decomposition of this system, can potentially 
be numerically inaccurate. The algorithm used 
is taken from [ Golub and Reinsch , (1971) ]. 

The orthogonality of the eigenvectors has been 
verified. However, with any linear vector space 
of this dimension, one must remain aware of 
possible numerical errors [Hamming, (1973) ; 

Wilkinson and Reinsch, (1971)]. For example, 
if the position of P or Q in an observation is the 
same as P or Q for another observation, the 
Gram matrix will not have independent 
elements. Though this is unlikely to ever to be 
exactly true, we have found that we must 
require all points P and Q be separated by 
more than a minimum distance. For large a 
system, n>1000, using a minimum distance of 
20 km is satisfactory. For systems with n = 
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several hundred, a minimum distance of 5 km works well. Figures Vn_2_f,g,h,i j, and 1 show a 
selection of the basis functions, x ¥ l , (6. 11) computed from these data and eigenvectors. The 
eigenvectors and basis functions have been ordered with the decreasing values of the 
eigenvalues: As expected from the general theory [ Parker , (1994) ] the x ¥ i 


n£ VII_2_b Spectrum 


Fig_VH_2_c Spectrum 


h=300 km 
separation = 3 deg 
log(detenninant)=-7035.0527 
I am (m ax)/ lam (min ) =3. 12el5 


h =450 Ion 
separation = 1 deg 
Iog(detcrmitum)=-7380.6212 
iam(max)/ iam(min )=4,71el3 


0 W W3 -KJ 


Fig_VH_2_d Spectrum 


m VII 2 e Grace Mission 


fa =450 km 
separation = 3 deg 
log(detenninai)t)=-7380.6212 
lam(max)/iam(min)=1.94el5 


GRACE Mission (1135 data) 
h=300 km 
separation = 3 deg 
log(deter!mjiant) =: -6927. 1294 
lam(max)/ lam(nrin)=7.57el3 


200 


300 


*» 


1000 


1200 














become more variable, oscillatory, as X, decreases in size. This supports the notion that by 
truncating the series of basis functions, , a certain smoothing is achieved by eliminating high 
frequency variations. 
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The selection of random positions and orientations, would be expected to provide the 
best set of n data points. Since the target of this analysis is the GRACE mission, the physical 
orientation of the data in low and mid latitudes will be controlled by the satellite inclination 
(87°) - near the pole, of course, all orientations will be obtained. Therefore, this experiment is 
repeated with the parameters given in Table VI 12 2. 


Table VI1_2_2 Parameters for GRACE Simulation 

Area Size 

20°x20° 

Geoid anomaly Size 

5°x5° 

Geoid anomaly amplitude 

8 cm 

Satellite Altitude 

400 km 

Satellite Separation 

3° (330 km) 

Number of Data Points 

1000 

Distribution of positions 

Random 

Distribution of orientation 

87°, 183° 


Figure VII_2_e shows the eigenvalues for this GRACE system, often referred to as the spectrum 
of the system. Table VII23 summarizes the eigenvalues of the cases presented. With the 
GRACE data we see the same range, 10 14 , of the eigenvalues is present as existed in the 


Table VII_2_3 Eigenvalue Summary 

h(km) 

Separation 

o 

Determinant 


^max 

300 

1.0 

5.05*1 O' 7076 

2.60*10 14 

3.24* 10 2 

300 

3.0 

8.86*1 O' 7036 

3.12*10 15 

2.24* 10 3 

450 

1.0 

2.39* 1 O' 7381 

4.71 * 10 13 

9.09* 10 1 

450 

3.0 

3.59* 10' 7334 

o 

X 

ON 

6.76* 10 2 

450 GRACE 

3.0 



1.00* 10 3 


previous case. Figures VII_2_m,n,o,p, and q show a selection of the basis functions, T, , (6 1 1 ) 
computed from these data and eigenvectors The same general character of the basis functions is 
evident The geometry of the GRACE data set does not seem to be a factor in application of GIT 
to geoid recovery. 
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VII3 Error analysis 

Errors in the measurement of potential difference come from two difference sources. 

First are the inherent measurement errors that will, for this analysis, be assumed to be 
uncorrelated, zero mean, Gaussian noise with a given standard deviation, a, i.e. normally 
distributed errors. The second source is the process of converting the observed velocity 
difference (A dp/dt) into the potential difference (AT). This involves both geometrical 
corrections (Section V) , correction for the dynamical interaction between AT and GM/r (section 
VII_5_3 below) , and other unmodeled orbital effects. The first source will be studied in this 
section. The second will be studied in following sections, and will provide some information 
about scaling the standard deviation to take these other phenomena in to account in the 
estimation process. 

The data sample described in Table VII_2_1 is used. Potential at satellite positions P,Q 
is computed using Poissons Integral Formula, (3.1) and the difference used as the observation. 
Gaussian noise is added, computed using an algorithm from [Hamming, (1973) , p 143] based on 
the system random number generator. These data are used to determine the constant C, (6. 1 7) 
based on the condition (6. 18) using the method of Regula Falsi [Hamming, (1973) , p 65], 

Recall, that the GIT SEM can obtain a solution for the expansion coefficients, that exactly 
matches the observations and has a minimum variance of the model. What we seek is the 
solution, that only matches the observations to within the observation uncertainty, and has the 
minimum model variance. For illustrative purposes, three observation uncertainties are used, as 
given in Table VTI_3_1 . Also given are the determined value of Cfk^, and C. Choosing C has 
the same practical effect as removing all \<~ C. 


h(km) 

Separation 

o 

n 

o (m/sec) 

C 

cfK~ 

300 

1.0 

1000 

1.0x1 O' 5 

6.44x10' 

1.98x1 O' 3 

300 

1.0 

1000 

1.0x10* 

6.73xl0‘ 2 

2.08x10* 

300 

1.0 

1000 

1.0x1 O' 7 

2.56x1 0' 3 

7.90x10* 

300 

3.0 

1000 

1.0x1 0' 5 

2.71x10' 

1.21x10* 

300 

3.0 

1000 

1.0x10* 

1.81xl0- 2 

8.08x10* 

300 

3.0 

1000 

l.OxlO' 7 

3.86x10* 

1.72x1 O' 7 


Figures VII_3_a, VII_3_b, and Vll_3 c, show the recovered geoid for these data sets. Clearly, 
higher accuracy, smaller a, provides greater resolution. 
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One important product of the GITSEM process is computation of a formal uncertainty of 
the model, (equation 6.19). To illustrate this, note that the model has zero signal outside the 5°x 
5° geoid anomaly. Figure VII_3_d, is the recovered geoid, where, a zero is plotted when the 
computed model is less than 3 times the computed uncertainty of the model. The values 
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predicted by the model outside the 5°x 5° geoid anomaly are not significantly different from 
zero, as determined by the model. 

Finally, the resolving power of the GITSEM, i.e. the sensitivity and resolution, is 
investigated as follows. With the analogue of a matched filter, we ask, what fraction of the 
signal amplitude is recovered? This will depend on the satellite altitude, the satellite separation, 
the size of the geoid anomaly, and the observation uncertainty. A series of simulations was done 
with the variables given in Table VI_3_2. 


Table VII32 Resolving Power of GITSEM 

Satellite Height (km) 

150 

200 

250 

300 

350 

400 

Satellite Spacing (km) 

110 

220 

330 



nn 

Anomaly Size (°) 

1 

2 

3 

4 

5 

Hi 

Data Accuracy (m/sec) 

1.0x1 O' 5 

l.OxlO -6 

1.0x1 O' 7 



■n 


Figures VII_3_e and VII_3_f summarize these results. Clearly lower satellite altitude and larger 


Figure VII 3 e Attenuation <3 dag spacing) 



: 1 60 tan ' 
700 km 
260 km , 
300 km 
350 km 
450 km | 


anomaly block size give higher resolution. As seen in numerous previous studies, there is an 
inherent limitation in recovery of geoid anomalies smaller than the satellite height. There is a 


- 45 - 



























smaller effect of satellite separation, larger separation providing somewhat more sensitivity. 


Figure VII_3_f Attenuation (1 <teg spacing) 



—•—1 SOtorTl 
—a— 200 Km | 
-a- 250 km 
-N- 300 km 
— •— 360 km 
—♦—450 km 
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VU_4 Simulation Geoid Block Anomalies 

An example is now given with several geoid block anomalies. Four anomalies are created 
as shown in the cartoon. Fig VTI 4 a and described in Table VII_4_1 . 


Figure VII_4_a 



Geometry of SST Solution 
32 cm Geoid height, 1x5 deg anomaly (twice) 
-12 cm Geoid height, 1 *10 deg anomaly (twice) 


Table VII 41 Block Anomaly Simulation 

height (km) 

300 

separation (°) 

1.0 

n (data) 

1000 

ct (m/sec) 

© 

X 

o 

Anomaly 1 (size,amplitude) 

1 ° x 5° ,32 cm 

Anomaly 2 (size,amplitude) 

1 0 *5° ,32 cm 

Anomaly 3 (size,amplitude) 

1°*10° ,-12 cm 

Anomaly 4 (size,amplitude) 

r*10° ,-12 cm 
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As before, the potential at satellite positions (P,Q) is computed using Poissons Integral Formula 
(3.1) and the difference used as the observation Gaussian noise is added The two 32 cm 
anomalies are resolved. However, only 75% of the amplitude (24 cm) is recovered. Similarly, 
the -12 cm anomalies are resolved with reduced amplitude (-8 cm). 
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VII 5 Orbit Simulation 


VI 151 Introduction 

The end to end orbit simulation process is described in figure VII51 . It begins with 
simulated satellite to satellite range rate data. The parameters of the simulations are summarized 
in Table VII_5_ 1 , and are intended to represent the GRACE mission. The initial state vectors are 
defined in Table 4 .1. However, before reviewing the simulation results, the calculation of the 
potential difference, AT, from the observed satellite to satellite range rate, dp/dt, must be 
discussed. Now, we are interested in the contribution of the anomalous potential, T, to the range 
rate. Therefore, the first step is to perform a Precision Orbit Determination (POD) using the SST 
data and any other available tracking data - here assumed to be ground based laser tracking - 
using a chosen reference potential, U, to obtain the orbit residuals, A dp/dt. These residuals 
contain the desired signal, the interaction terms from n5r, and other orbital errors due to 
inadequate force models. Therefore, we now discuss the following three topics: 1) The filter 
parameters, 2) Calculation of the interaction correction, nSr, and 3) the geoid database used for 
calculation and comparison of the results. 



Table V3I 51 Simulation Parameters 

Data Geometry 

400 km altitude, 1=87°, 400 km separation 

SST Data 

10 second range rate data 
60 Day mission 

Ground Tracking Stations 

1 8 laser ranging stations tracking both satellites 

Force Model 

Total Geopotential: 1=2,70, 2-90, 2-180 
Sun, Moon Tides (Ocean Body) 

Drag (MSIS), Solar Pressure, Earth Albedo 
Reference Potential: 1=2-40, 2-90 

Orbit Arcs 

1 day 

Solve for parameters: drag scale factor 

Map Region 

Eastern Mediterranean 

Filter parameters 

X, ow =800 sec, ^=30 sec, (5600 km, 47 °) 


VII 5 4 Filter Parameters 

To calculate the potential difference (AT) from the observed satellite to satellite tracking 
velocity residuals (A dp/dt) we start with equation (4.33) is used. Recall, the AT response of the 
satellite velocity was local, and takes place immediately. However, the interaction perturbation 
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Data Processing Structure 


Tracking & SST Data 


Filter bounds 


4>i. K 


o i , t 
J ram 


Precision Orbit 
Determination 




/ 


n8r has an instantaneous component, but also has a cumulative effect, that results in periodic 
changes with, in fact, the orbital period Since we are interested in the instantaneous change due 
only to AT we chose to filter out signals with periods longer than the time passing over the 
region of interest. The satellite traveling at 7.5 km/sec, passes over 40° of latitude (4400 km) in 
less than 600 seconds. So filtering signals with periods longer than 600 seconds or more, would 
leave all the instantaneous signals. In fact, the results are not very sensitive to the period cutoff, 
we experimented with cutoff 5 s of 600, 800, 1200, 1600, 2400, and 4800 seconds, and adopted 
800 seconds for all the analysis. The filtering was accomplished using a Fast Fourier Transform 
Algorithm [Press et al., (1991) ]. Each one day arc of residuals (8640 data), was filtered 
separately. The data were padded with zeros to obtain a sequence of exactly 2" points. The 
Fourier coefficients were then convolved with a 10 lh order Butterworth filter. To mitigate 
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measurement noise, i.e. high frequency noise, a 30 second short period cutoff was also applied. 
The resulting Fourier coefficients were then converted back to the time domain. 

Figures VII_5_2_a and VII_5_2_b illustrate the result of this filter. The test case has a total 



geopotential (U+T) to degree and order 70, and a reference geopotential (U) to degree and order 
40. Figure VII 5_2_a shows the computed residual (A dp/dt) and the filtered values (A dp/dt sr „ u , n ) 
for three passes over the chosen region. Figure VII_5_2_b shows the filtered values of the 
computed residual (A dp/dt**^) , the computed values of the potential difference (AT), and the 
difference, and the difference (A dp/dt^.^ -AT). For this comparison, the computed potential 
difference is converted to velocity using 4 > =TJ+T=yN= r v8v. The variances of the computed 
residuals (0.0252 cm/sec = 252 p/sec) and the computed potential difference (0.0246 cm/sec = 
246 p/sec) are nearly the same. The standard deviation of the difference (0.0024 cm/sec) is about 
10% of the signal, as predicted in section III. The difference between the computed residual and 
the potential difference is due to the interaction term, to be discussed in the next section, and 
other orbital errors. 
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VII 5 3 Calculation of n5r 


If one knows the anomalous potential (T), then we have shown that the correction (n8r) 
can be computed with sufficient accuracy. Here, this is not the case. Therefore, we develop an 
approximation based on the available data. Recall that the computed residuals A dp/dt 
~AT/v=(T(P )-T(Q))/v. Therefore, if the two satellites are in the same circular coplanar orbit then 
one can use the computed residuals, A dp/dt, to compute an approximation to the potential along 
the orbit. This is only approximately true because of the inevitably different evolution of the 
orbits, and that the earth rotation will result is a different potential at the same orbit position at a 
different time. We give in Figure VII_5_3_a the potential recovered by this method, expressed as 
velocity, and the computed anomalous potential, T. Again the case studied is total potential 
T=U+T complete for 1=2-70, the reference potential U:2-40, and the anomalous potential T:41- 
70. Now, the correction 5r depends on the combination 2T+r(dT/<9r). We obtain the second term 
as follows. Assume that we have the spherical harmonic expansion for T, i.e. 


T(r,(p.l) 


GM 

r 


I- 1 max 

I 

l- / min 



m~ l _ 

X ^L(sin <p)(C lm cos ml + S lm sin ml) 

m~ 0 


(7.1) 
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(7.2) 


T(r, <p, A) = ]T T, (r, <p, X ) 

( T 

and the desiderata 


JT 

r — 


I ( / + 1 )4tT, 

/=/min ' 


(7.3) 


Consider, an arc, in any orientation on the earths surface containing the point of interest (<p 0 ,X. 0 ) in 
the orbit subsatellite track. The arc, and (tp 0 ,X. 0 ) define a coordinate system and the spherical 
harmonics can be transformed to this system. Now, along this arc, the longitude is zero and (7.1) 
reduces to an expansion in Legendre polynomials. So we take the recovered T, along the arc, and 
expand it in Legendre polynomials through degree 142. This expansion is centered on the point 
of interest, ((p 0 ,/. 0 ) , and for -jt/2 to n/2 along the orbit. From this expansion, the sum (7.3) is 
evaluated at (<p„,L 0 ) Figure VII 5 3 b shows this gradient is shown as recovered from A dp/dt, 
and computed from the, in this case, known potential. The same tracks are used as before. 
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The generally good agreeement between the recovered and model values of r(c?T/dr) is 
satisfactory considering the purpose here is to calculate a 10% correction. The correction is 
calculated for each satellite separately, and then n8r=n(8r(P)-8r(Q)) is computed and added to the 
observed A dp/dt to obtain the estimate of AT/y to be analyzed. Recall, that the calculated Sr will 
introduce long period perturbations. Therefore, the result for 8r, will also be filtered before being 
used, using a 6000 second cutoff period. Figure VII_5_3_c shows the computed values for n8r(P) 
and of n8r. nSr(Q) is not plotted simply because it is shifted in time from n8r(P) by p/v ~ 400 
(km)/ (7.5 km/sec) = 53 seconds. 
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time (1 Q*aec) 


How well does this calculation work? In figure Fig_VII_53_d we show the filtered orbit 
residual Adp/dt minus the computed dv(model)/dt, and the calculated n8r, for the same three 
tracks used above. One sees the signature of ndr. In addition there is a larger signal due, 
presumably to other unmodeled effects. Figure VIl_5_3_e shows the result of correcting the 
observed velocity with n8r. 

Examination of many samples of data, give similar results. It is clear, that the conversion 
of observed range rate residual, Adp/dt, to AT/y has remaining unmodeled errors. Though only of 
the order of 5% to 10% of the desired signal, they limit the accuracy of the recovered geoid, as 
will be discussed below. 
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VI_5_5 Geoid Data Base 

For facilitating comparison of the recovered geoid, N, with the model geoid, N(model), a 
geoid data base has been created. The fundamental coordinates used are geodetic latitude and 
longitude. The surface of the earth is divided into 0.5° (55 km), approximately equal area blocks. 
These blocks are arranged in 0.5° latitude bands, with coordinate <p=0.25 o ±kx0.5°. Each latitude 
band is divided into n=[cos(<j> 0 )x 360.0/0.5 + 0.5] blocks, where the symbol [x] signifies the 
integral part of x. The longitude blocks have the size A>.=360.0°/n. When extracting data in the 
data base, there are two options. First, is to return the same value for all points within the block. 
Second, is to use linear interpolation, as a ruled surface, for a point bounded by four data points. 

When creating a geoid data base from the potential expressed in spherical harmonics, e g. 
the EGM96, we assume that the independent variables are geocentric latitude and longitude. 
When evaluating the spherical harmonics, for insertion in the geoid data base for a geodetic 
latitude and longitude, these are converted to geocentric latitude and longitude for the 
calculation. Note that this is consistent with generation of satellite observation residuals, Adp/dt. 
In archiving the satellite residuals the coordinates of the two satellites (P,Q) are archived as 
height above the reference ellipsoid, geodetic latitude and longitude, as well and the geocentric 
Cartesian coordinates and velocities. 
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VII 6 Orbit Simulation 


VII 61 Introduction 

In previous sections, necessary elements of. Potential Theory (III), Orbit Theory (IV), and 
Geophysical Inverse Theory ( VT> to analyze satellite to satellite tracking data for geoid 
determination have been reviewed. Some numerical experiments were described that illustrate 
how these elements can be combined, and what sort of results may be expected. These numerical 
experiments were progressively more complex and complete First the nature of the eigenvalues 
and eigenfunctions that come from the singular value decomposition were described. These 
results depend only on the physical geometiy of the measurement, and is independent of the 
measurement itself. Second, some properties of the solution assuming a direct measurement of 
potential difference: accuracy, resolution, error propagation, and sensitivity were shown. For 
these analyses, idealized geoid anomalies were used: blocks with a dimension of 1 .0 degree and 
5.0 degrees were combined. Now a number of orbital simulations will be offered, with increasing 
complexity in the desired geoid, and other orbital effects. For these the geoid model will be the 
EGM96 geopotential model [Lemome and al., (1998) ]. EGM96 is a complete Spherical 
Harmonic representation to degree and order 360. For the illustrative proposes of these 
demonstrations, the EGM96 model will be truncated, at degree 70, 90 and 180. In addition other 
physical forces will be included: lunar and solar forces, ocean and solid body tides, atmospheric 
drag, solar radiation pressure, and earth albedo pressure. The focus of these last simulations will 
be to demonstrate the recovery of geoid information from satellite to satellite range rate, which 
involves correcting for the indirect effects (equation 4.33) in the presence of other force model 
errors, that can only approximately be accounted for in the orbit fit. Exploration of the 
degradation of the recovered geoid in the presence of these other forces, is a much larger effort 
than was addressed here, but some limited cases will be shown. 

The parameters of the simulations are summarized in Table VII51, and the initial state 
vectors are defined in Table 4. 1 . These are intended to represent the GRACE Mission. The 
end-to-end orbit simulation process is described in Figure VII_5_1 . It begins with simulated 
satellite-to-satellite tracking (SST) range rate data. Now, we need to obtain the contribution of 
the anomalous potential, T, from the observed range rate. Therefore, the first step is to perform a 
Precision Orbit Determination (POD) using the SST data and any other available tracking data - 
here assumed to be ground based laser tracking - using a chosen reference potential, U, to obtain 
the orbit residuals, A dp/dt. These residuals contain the desired signal, the interaction terms from 
n6r, and other orbital errors due to inadequate force models. Therefore, the calculation of the 
potential difference, AT, from the observed satellite to satellite range rate, dp/dt, is accomplished, 
as described above, using the POD residuals. A number of cases are computed, summarized in 
Table VII 6 1. 
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Table VII 61 Parameters of Orbit Simulations 

Case 

Total 

Potential 

U+T 

Reference 

Potential 

U 

Other Forces 

Latitude 

Range 

Longitude 

Range 

Sample 
POD a 
(cm/sec) 

D 

1=2-70 

1=2-40 

Sun,Moon,Tides 

20°-50° 

20° -30° 

0.03 

2 

1=2-180 

1=2-90 

Sun,Moon,Tides 

20°-50° 

20° -30° 

0.01 

3 

1=2-90 

1=2-40 

Sun,Moon,Tides 
Drag, Solar 
Pressure, Earth 
Albedo 

0 
0 

1 

C/i 

o 

0 

20° -30° 

0.03 

4 

1=2-180 

1=2-90 

Sun,Moon,Tides 
Drag, Solar 
Pressure, Earth 
Albedo 

20° -50° 

20°-30° 

0.01 

5 

1=2-180 

1=2-90 

1 

Sun,Moon,Tides 
Drag, Solar 
Pressure, Earth 
Albedo 

29°-41° 

230-270 

0.01 


VI16 2 Region of Interest. 

The main thrust of using the GITSEM is to obtain local or regional maps of the geoid, 
from the measurement of SST range rate We have chosen one region for the simulation: the 
Eastern Mediterranean Sea. This region is tectonically active, and has a significant signal at all 
wavelengths. Therefore, it provides a useful geophysical test bed for the method. Figure VIl_6_a 
is a free air gravity anomaly map of the region, centered roughly on the island of Crete. 

VII 6 3 Case 1 

We begin with the total potential, U+T, complete for degree and order 1=2 through 70 and 
the reference potential, U, complete for 1=2 though 40. Therefore, the anomalous potential 
sought, is complete for 1-41 though 70. The 60 day mission is processed, and the observations 
selected requiring that no positions P or Q are closer than 20 km. This resulted in 1187 
observations. The Gram Matrix calculation and Singular Value Decomposition provided the 
eigenvalues and eigenvectors. Smoothing was achieved by selection of the Lagrange Multiplier, 
C. Table VII 6 2 gives the values. Figure VLI_6_b shows the recovered geoid, and Figure 
VTI_6_c shows the anomalous geoid (1=41 through 70). Also given in Table VII 6 2 are 


- 59 - 






































• (d*g) 


Figure VII_6_a 


Eastern Mediterranean 



- 60 - 






- 61 - 











statistics of the geoid fit. We see the general agreement of features. The statistics of this solution 
is as follows. The variance of the geoid for the selected region is 268.0 cm. The rms of the geoid 
fit is 88.2 cm. To identify the error sources in this map, the orbit positions (P,Q) are used, i.e. the 
same Gram Matrix, and the observed potential difference is calculated from the geoid model 
using the Poisson Integral formula (3. 1 ). For this calculation the rms of the geoid fit is 23.0 cm 
Figure VII_6_d shows this recovered geoid. This calculation is repeated selecting orbit positions 
with a minimum position difference of 10 km. This results in n=l 881 observations, and has an 
rms of geoid fit of 16.7 cm. This is a marginal improvement considering the increased dimension 
of the basis functions. A number of remarks are in order. First, it is not possible to obtain exact 
agreement between the spherical harmonic representation and the G1TSEM. This should not be 
surprising, considering the fact that these use fundamentally different basis functions, and they 
are finite in dimension. Second, the GITSEM does have the potential to achieve centimeter 
accuracy in geoid height. Third, errors in the recovered geoid are due to errors in converting the 
observed range rate residuals, A dp/dt, to potential difference,AT. 


- 62 - 



Table VTI_6_2 Summary of Geoid Results 


Case 

n 

^Srnim^^^tnax 

Obs a 
(cm/ sec) 

Var(Obs) 
(cm/ sec) 

Geoid 
height a 
(cm) 

Var 

(geoid) 

(cm) 

Comment 

■ 

1187 

c/^io - 4 

8.9 1 x 1 0' 1 

21.85 

88.2 

268.0 

20 km 
separation 

■ 

1187 

CA^-10' 7 

3.61x1 0' 3 

21.97 

23.0 

268.0 

Upward 

Continuation 

■ 

1881 

C/^=10- 9 

3.18xl0 4 

20.90 

16.7 

268.0 

10 km 
separation 

B 

1186 

1 0' 2 ( 1 07) 

7.59* 10- 2 

2.67x10-' 

64.37 

86.38 

20 km 
separation 

2 

1186 

10" 6 (402) 

1 .61 * 10 -4 

2.54x10-' 

21.17 

86.38 

Upward 

Continuation 

2 

1913 

1 O' 2 ( 1 3 1 ) 

7.90*1 O’ 2 

2.68x10' 

69.72 

86.38 

10 km 
separation 

3 



8.90x1 O' 1 

22.59 

85.52 

303.2 

20 km 
separation 

B 



2.92 x 10" 2 



86.38 

20 km 
separation 

5 

271 

C/X^-IO* 

3.01 xio- 2 

3.36x10-' 

53.1 

124.5 

5 km 

separation 
Limited Area 

VII 6 4 Case 2 


In this case we illustrate recovery of smaller scale features of the geoid. Therefore, we 
begin with the total potential, U+T, complete for degree and order 1=2 through 180 and the 
reference potential, U, complete for 1=2 though 90. Therefore, the anomalous potential sought, is 
complete for 1-91 though 180. As before, the 60 day mission is processed, and the observations 
selected requiring that no positions P or 0 are closer than 20 km. This resulted in 1186 
observations. The Gram Matrix calculation and Singular Value Decomposition provided the 
eigenvalues and eigenvectors. Smoothing was achieved by selection of a minimum eigenvalue. 
Table VII62 gives the values. For the cases where X^ has been chosen, the number of 
eigenfunctions used for the model is given in parenthesis. Figure VII_6_e shows the recovered 
geoid, and Figure VII_6_f shows the anomalous geoid (1=91 through 180). Also given in Table 
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VII62 are statistics of the geoid fit. We see the general agreement of features. The statistics of 
this solution are as follows. The variance of the geoid for the selected region is 86.38 cm The 
rms of the geoid fit is 64.37 cm. To identify the error sources in this map, the orbit positions 
(P,Q) are used, i.e. the same Gram Matrix, and the observed potential difference is calculated 
from the geoid model using the Poisson Integral formula (3. 1 ). For this calculation the rms of the 
geoid fit is 2 1 . 1 7 cm. Figure V7I_6_g shows this recovered geoid. 



VII 6 5 Case 3 

In this case we illustrate geoid recovery in the presence of non gravitational forces, 
thermosphere drag, solar radiation pressure, and earth albedo pressure. In this case the MSIS 
thermosphere model is used, a solar pressure model, and an earth albedo pressure model, in 
generating the simulated data, and in the precision orbit determination (POD). In the POD, a 5% 
error is introduced in the overall model, and a drag scale “solve for” factor is included. We 
begin with the total potential, U+T, complete for degree and order 1=2 through 90 and the 
reference potential, U, complete for 1=2 though 40. Therefore, the anomalous potential sought, is 
complete for 1-41 though 90. As before, the 60 day mission is processed, and the observations 
selected requiring that no positions P or Q are closer than 20 km. This resulted in 1135 
observations. The Gram Matrix calculation and Singular Value Decomposition provided the 
eigenvalues and eigenvectors. Smoothing was achieved by selection of a minimum eigenvalue. 
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a,,,,,,. Table VII 6 2 gives the values. For the cases where has been chosen, the number of 
eigenfunctions used for the model is given in parenthesis. Figure VII 6 h shows the recovered 
geoid, and Figure VTI_6_i shows the anomalous geoid (1=41 through 90). Also given in Table 
VTI62 are statistics of the geoid fit. We see the general agreement of features. The statistics of 
this solution as as follows. The variance of the geoid for the selected region is 303.2 cm. The 
rms of the geoid fit is 85.52 cm. The fact that the results here are so similar to Case 1 , is 
encouraging, in that the presence of nongravitational forces, and errors in modelling them, seem 
to have a small effect on the geoid recovery. This is consistent with the idea that the effects of, 
and errors in, drag and solar pressure have quite a different spectral character, and should be 
separated from the geoid. However, the thermosphere models used here are quite smooth, and do 
not even attempt to represent small scale thermosphere variations such as gravity waves, winds, 
etc. No doubt, great care will be necessary in treating these non-gravitational forces. 

Vll_6_6 Case 4 

In this case we continue to illustrate geoid recovery in the presence of non gravitational 
forces, thermosphere drag, solar radiation pressure, and earth albedo pressure, but for finer 
scale geoid features. The MSIS thermosphere model, a solar pressure model, and an earth 



albedo pressure model, are used generating the simulated data, and in the precision orbit 
determination (POD). In the POD, a 5% error is introduced in the overall model, and a drag 
scale “solve for” factor is included. We begin with the total potential, U+T, complete for 
degree and order 1=2 through 1 80 and the reference potential, U, complete for 1=2 though 90. 
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Therefore, the anomalous potential sought, is complete for 1-91 though 180. As before, the 60 
day mission is processed, and the observations selected requiring that no positions P or Q are 
closer than 20 km. This resulted in 1132 observations. The Gram Matrix calculation and 
Singular Value Decomposition provided the eigenvalues and eigenvectors. Smoothing was 
achieved by selection of a minimum eigenvalue, Table VII_6_2 gives the values. For the 
cases where has been chosen, the number of eigenfunctions used for the model is given in 
parenthesis. Figure VII_6J shows the recovered geoid. Also given in Table VII_6_2 are 
statistics of the geoid fit. We see the general agreement of features. The statistics of this 
solution as as follows. The variance of the geoid for the selected region is 86.38 cm. The rms 
of the geoid fit is 57.01 cm. 



VII 6 7 Case 5 


In this case we continue to illustrate geoid recovery in the presence of non gravitational 
forces, thermosphere drag, solar radiation pressure, and earth albedo pressure, for finer scale 
geoid features. Here, we reduce the area of interest, and decrease the minimum distance 
between points (P or Q), to see if this would increase the recovered geoid resolution. As above, 
the MSIS thermosphere model, a solar pressure model, and an earth albedo pressure model, are 
used generating the simulated data, and in the precision orbit determination (POD). In the 
POD, a 5% error is introduced in the overall model, and a drag scale “solve for” factor is 
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included. We begin with the total potential, U+T, complete for degree and order 1=2 through 
1 80 and the reference potential, U, complete for 1=2 though 90. Therefore, the anomalous 
potential sought, is complete for 1-91 though 180. As before, the 60 day mission is processed, 
and the observations selected requiring that no positions P or Q are closer than 5 km. This 
resulted in 271 observations. The Gram Matrix calculation and Singular Value Decomposition 
provided the eigenvalues and eigenvectors. Smoothing was achieved by selection of a 
minimum eigenvalue, A^. Table VII62 gives the values. For the cases where k mm has been 
chosen, the number of eigenfunctions used for the model is given in parenthesis. Figure 
VII 6_k shows the recovered geoid and Figure VII_6_1 shows the anomalous geoid (1=91 
through 180). Also given in Table VII62 are statistics of the geoid fit. We see the general 
agreement of features. The statistics of this solution as as follows. The variance of the geoid 
for the selected region is 124.5 cm. The rms of the geoid fit is 53. 1 cm. 
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VIII Discussion 


The geoid recovered in all cases represents the major features of the model geoid. 

It is not possible to obtain exact agreement between the spherical harmonic representation and 
the GITSEM. This should not be surprising, considering the fact that these use fundamentally 
different basis functions, and they are finite in dimension. Second, the GITSEM does have the 
potential to achieve centimeter accuracy in geoid height. Third, errors in the recovered geoid 
are due to errors in converting the observed range rate residuals, A dp/dt, to potential 
difference, AT 

The present accuracy of the orbital simulation, does not meet the objective of achieving 
0.5 p/sec range rate residuals, and centimeter accuracy short wavelength geoid height.. Due to 
limited scope of this analysis, refinements of calculating the potential difference AT, from the 
range rate residuals, A dp/dt, could not be explored. The current method has an error of ~5%. 
There are a number of avenues to be investigated. First would be use of short arc calculation to 
obtain the residuals. Since we believe that the remaining error is due to unmodeled orbital 
errors, this could reduce these errors. The second would be to use an iterative method. Taking 
the first solution for the geoid, calculate a correction to the residuals, and repeat the GIT geoid 
estimate to obtain a correction. A third method would be to develop a recursive digital filter, 
based on the orbit model, to improve the correction estimate. 

In a careful and detailed analysis Jekeli, [Jekeli, (1981) ] has shown two critical results. 
Jekeli combines the error of representation (3.20) and the error of omission (3.19 suitably 
generalized for r<R as 3. 17) as the downward continuation error - the error of commission is 
not included. Jekeli, with « = 300 , finds (tables 3 and 4) that the DCE largest near the poles 
o(e r ) <0.090 mgals (0.290 mgals max) gravity anomaly and o(e r ) <0.042 cm (0. 14 cm max) 
geoid height. Second, by seeking anomalies averaged over a spherical cap of about 1 .4 degrees 
(tables 6 and 7), the DCE - again largest at the poles, is estimated to be o(e r ) <0.004 mgals 
(0.014 mgals max) gravity anomaly and c(e r ) <0.0020 cm (0.0066 cm max) geoid height. 
Therefore the conclusion of Jekeli [Jekeli, ( 1981) , p 127] “The downward continuation errors 
depicted in tables 3 through 7 are completely insignificant with respect to anticipated 
measurement accuracies of 1 mgal and 10 cm in the gravity anomaly and geoid undulation, 
respectively.” And, “ .. the estimation of point or mean gravity anomalies and geoid 
undulations (height anomalies) using the outer series expansion to degree 300 anywhere on the 
earth’s surface is practically unaffected by the divergence of the total series.” We believe that 
using a spherical approximation for the anomalous potential, has sufficient accuracy for our 
geoid estimate. 

We turn now to the question of computer resources needed for this analysis. This study 
has been done on an INTEL Pentium II, 233 Mhz PC. Storage on this class of machine is quite 
adequate. However, the processor speed limits the analysis. There are two major steps in the 
computation: 1 ) the Precision Orbit Determination (POD) and 2) calculation of the Gram 
Matrix. For this analysis, we use the canonical 60 day mission. The POD is accomplished in 
less that 1 day, the exact amount of time, of course, depending on the force model complexity. 
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and the number of iterations (i.e. on the amount of bad data to be screened). This is not 
limiting. The algorithms developed for calculating the Gram Matrix, require about 0.2 sec per 
matrix element on this machine. Being a symmetric matrix the time for each region would be 
0.2 xnx(n+ 1 )/2 seconds. So for n= 1 000, we need about 28 hours of processor time. This is the 
bottleneck. Assuming the regional dimension used in this analysis, 20x30=600 square degrees, 
there would be 70 such regions for the whole earth. Presumably, for global analysis, one would 
want to use the classical spherical harmonic analysis anyway. In addition, the PC processor 
used here is now obsolete: there are 1 Ghz processors on the market now. A factor of 10 
improvement in running time is easily achieved with PC architecture today, and greater speed is 
available with more advanced and expensive computer architectures. If desired, once set up, 
one could produce a global set of geoid maps in a few days. It seems that this approach is 
feasible. 
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